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ABSTRACT 

We study the orbital parameters distribution of stars that are scattered into nearly radial orbits and then spiral 
into a massive black hole (MBH) due to dissipation, in particular by emission of gravitational waves (GW). This 
is important for GW detection, e.g. by the Laser Interferometer Space Antenna (LISA). Signal identification 
requires knowledge of the waveforms, which depend on the orbital parameters. We use analytical and Monte 
Carlo methods to analyze the interplay between GW dissipation and scattering in the presence of a mass sink 
during the transition from the initial scattering-dominated phase to the final dissipation-dominated phase of 
the inspiral. Our main results are (1) Stars typically enter the GW-emitting phase with high eccentricities. (2) 

The GW event rate per galaxy is few x 10 9 yr -1 for typical central stellar cusps, almost independently of the 
relaxation time or the MBH mass. (3) For intermediate mass black holes (IBHs) of ~ 10 3 M 0 such as may 
exist in dense stellar clusters, the orbits are very eccentric and the inspiral is rapid, so the sources are very 
short-lived. 

Subject headings: black hole physics — stellar dynamics — gravitational waves 


1. INTRODUCTION 

Dissipative interactions between stars and massive black 
holes (MBH s; M. > 10 6 M 0 ) in galactic nuclei (e.g. Geb- 
hardt et al. [200(1 l200l . or intermediate mass black holes 
(IBHs; 10 1 2 < M. < 10 4 M 0 ), which may exist in dense stellar 
clusters, have been in the focus of several recent studies. The 
interest is mainly motivated by the possibility of using the dis¬ 
sipated power to detect the BH, or to probe General Relativity. 
Examples of such processes are tidal heating (e.g. Alexander 
& Morris l2003t Hopman, Portegies Zwart & Alexander 2004) 
and gravitational wave (GW) radiation (Hils & Bender 1995; 
Sigurdsson & Rees ll 997t Ivanov 120021 Freitag l200Tl 120031) 

A statistical characterization of inspiral orbits is of inter¬ 
est in anticipation of GW observations by the Laser Inter¬ 
ferometer Space Antenna (LISA). LISA will be able to ob¬ 
serve GW from stars at cosmological distances during the 
final, highly relativistic phase of inspiral into a ~ 10 6 A/ 0 
MBH, thereby opening a new non-electromagnetic astronom¬ 
ical window. GW from inspiraling compact objects (COs) is 
one of the three major targets of the LISA mission (Barack & 
Cutler 12003 . 2004 ; Gair et al. 12004 ). together with cosmolog¬ 
ical MBH-MBH mergers and Galactic CO-CO mergers. 

LISA can detect GW emission from stars with orbital period 
shorter than Pl ~ 10 4 s. In order for the shortest possible 
period to be small enough to be detectable by LISA , the MBH 
has to be of moderate mass, M. < 5 x 10 6 M 0 (Sigurdsson & 
ReesdU. LISA is expected to be able to detect inspiral into 
MBHs of 10 6 Mq to distances as far as > 1 Gpc. 

The detailed time-evolution of the GW depends on the ec¬ 
centricity of the stellar orbit, and therefore probes both Gen¬ 
eral Relativity and the statistical predictions of stellar dynam¬ 
ics theory. Due to the low signal to noise ratio, knowledge 
of the wave forms is required in advance. For this purpose, 
it is necessary to estimate the orbital characteristics of the 
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GW-emitting stars, and in particular the distrib ution function 
(DF) of their eccentri cities ( Pierro et al. 1200 It Glampedakis, 
Hughes & Kennefick 120021) . as the wave forms are s trong 
functions of the eccentricity (e.g. Barack & Cutlcr f2003i : Wen 
& G air 120051 ). This study focuses on inspiral by GW emis¬ 
sion. However, it should be emphasized that inspiral is a gen¬ 
eral consequence of dissipation, and the formalism presented 
below can be extended in a straight-forward way to other dis¬ 
sipation processes, such as tidal heating. 

The prompt infall of a star in to a MBH and its destruction 
have been studied extensively ( T2.H . Here we analyze a dif¬ 
ferent process, the slow inspiral of stars (Alexander & Hop- 
man 12003). A star on a highly eccentric orbit with small pe- 
riapse r p , repeatedly loses some energy A E every periapse 
passage due to GW emission, and its orbit gradually decays. 
At a distance r o from the MBH, where the orbital period is / J ( j, 
the time-scale to for completing the inspiral (i.e. decaying to 
aP^O orbit) is much longer than the time-scale If needed 
to reach the MBH directly on a nearly radial orbit. While the 
orbit decays, two-body scatterings by other stars continually 
perturb it, changing its orbital angular momentum J by order 
unity on a timescale tj. Because to If, inspiraling stars 
are much more susceptible to scattering than those on infall 
orbits. If to > tj , either because ro is large or r v is large 
(small A E), then the orbit will not have time to decay and 
reach an observationally interesting short period. Before that 
can happen, the star will either be scattered to a wider orbit 
where energy dissipation is no longer efficient, or conversely, 
plunge into the MBH. Inspiral is thus much rarer than direct 
infall. The stellar consumption rate, and hence the properties 
of the stellar distribution function (DF) at low J, are domi¬ 
nated by prompt infall, with inspiral contributing only a small 
correction. This DF describes the parent population of the 
inspiraling stars. 

We show below that the DF of the small subset of stars on 
low-J orbits that complete the inspiral and are GW sources 
is very different from that of the parent population (Fig. |3- 
This results from the interplay between GW dissipation and 
scattering in the presence of a mass sink during the transition 
(to ~ tj) from the initial scattering-dominated phase to the 
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final dissipation-dominated phase of the inspiral. 

This paper is organized as follows. In El we recapitulate 
some of the results of loss cone theory for the prompt infall, 
and extend it to slow inspiral. In El and 0we present a de¬ 
tailed analytical discussion of the main effects that determine 
the rates of GW events and their statistical properties. In El 
we describe three different approaches for studying the prob¬ 
lem: by Monte Carlo simulations, by solving numerically the 
2D diffusion / dissipation equation in E and J, and by a sim¬ 
plified analytical model, which mimics the behavior of a typ¬ 
ical star. In El we apply the MC simulation to a MBH in a 
galactic nucleus, and an IBH in a stellar cluster. We summa¬ 
rize our results in El 


2. THE LOSS-CONE 

The rate at which stars are consumed by an MBH and 
the effect this has on the stellar DF near it have been stud- 



N-body simulations with stellar captures were recen tly per- 
formed by Baumgardt, Makino, & Ebisuzaki ( 2004 ai 2 004 h) 
and by Preto, Merritt & Spurzem (20Q4j). 

We begin by summarizing these results, neglecting dissi¬ 
pative processes. We then extend the formalism to include 
dissipative processes. 

2.1. Prompt inf all 

The stellar orbits are defined by a specific angular momen¬ 
tum J and relative specific energy e = ip(r) — v 2 /2 (hereafter 
“angular momentum” and “energy”), where ip is the relative 
gravitational potential, and v is the velocity of a star with re¬ 
spect to the MBH. A spherical mass distribution and a nearly 
spherical velocity distribution are assumed. 

Orbits in the Schwarzschild metric (unlike Keplerian orbits) 
can escape the MBH only if their angular momentum is high 
enough, J > Ji c (e ). The phase space volume J < Ji c is known 
as the “loss-cone”. As argued below, stars that are scattered 
to low-J orbits are typically on nearly zero-energy orbits. For 
such orbits 


t r ~ e/e, whereas diffusion in J-space occurs on the angu¬ 
lar momentum relaxation timescale, 

tj~ J 2 /(J 2 )~ [J/J m (e)ftr, (2) 

where J m (e) is the maximal (circular orbit) angular momen¬ 
tum for specific energy e.The square root dependence of J 
on tj reflects the random walk nature of the process. Typi¬ 
cally, Ji c <C J m - In principle, stars can enter the loss-cone, 
J < Jicif) either by a decrease in J, or by an increase in e 
(up to the last stable orbit). In practice, diffusion in J-space 
is much more efficient: the energy of a star must increase by 
many orders of magnitude in order for it to reach the loss- 
cone, which takes many relaxation times. The angular mo¬ 
mentum of the star, on the other hand, needs only to change 
by order unity in order for the star to be captured, which hap¬ 
pens on a much shorter time t j < t r . 

The ratio between J\ c and the mean change in angular mo¬ 
mentum per orbit, A J, defines two dynamical regimes of loss- 
cone re-population (Lightman & Shapiro fl 977l) . In the “Dif¬ 
fusive regime” of stars with large e (tight orbits), A J <C Ji c 
and so the stars slowly diffuse in J-space. The loss-cone 
remains nearly empty at all times since any star inside it is 
promptly swallowed. At J > the DF is nearly isotropic, 
but it falls logarithmically to zero at J > J; c (Eq. 1101 , In 
the “full loss-cone regime” (sometimes also called the “pin¬ 
hole” or “kick” regime) of stars with small e (wide orbits), 
A J Jic and so the stars can enter and exit the loss-cone 
many times before reaching periapse. As a result, the DF is 
nearly isotropic at all J > J/ c . We argue below that only the 
diffusive regime is relevant for inspiral. 

The DF in the diffusive regime is described by the Fokker- 
Planck equation. We follow Lightman & Shapiro J1977l) , who 
neglect the small contribution of energy diffusion to J 7 , and 
write the Fokker-Planck equation for the number density of 
stars N(e,J; t ) as 1 

dN(e,J;t) dF{e;J) 

at dj ’ 1; 

where 

J) = N(s, J)(AJ) - J)(AJ 2 ). (4) 

The diffusion coefficients (A J) and (A J 2 ) obey the relation 


< 4:0 = 0) 


4 GM, 
c 


( 1 ) 


The size of the loss-cone J; c is nearly constant over the rel¬ 
evant range of e. Only during the very last in-spiral phase 
the energy of the star becomes non-negligible compared to 
its rest-mass, in w hich case the loss-cone is slightly modified 
(see section 14.11 1. Deviations from geodetic motion due to 
tidal interactions are neglected here. This assumption is justi¬ 
fied for COs orbiting ~ 10 6 M 0 MBHs, where the tidal radius 
is much smaller than the event horizon. For main-sequence 
(MS) stars where the tidal radius lies outside the event hori¬ 
zon, the loss-cone is similarly defined as the minimal J re¬ 
quired to avoid tidal disruption. 

Stars that are initially on orbits with J < will promptly 
fall into the MBH on an orbital timescale. Subsequently, the 
infall flow in J-space, E(e;J), is set by the rate at which 
relaxation processes (here assumed to be multiple two-body 
scattering events) re-populate the loss-cone orbits. 

Diffusion in e-space occurs on the relaxation timescale. 


(A J 2 ) = 2 J (A J) (J«J m ), (5) 

(Lightman & Shapiro [197 7: Magorrian & Tremaine 11999 1. 
Much of the difficulty in obtaining an exact solution for the 
Fokker-Planck equation stems from the dependence of the dif¬ 
fusion coefficients on the DF; self-consistency requires solv¬ 
ing a set of coupled equations. For many practical applica¬ 
tions the diffusion coefficients are estimated in a non-self- 
consistent way, for example by assuming local homogeneity 
and isotropy (e.g. Binney & Tremaine l 19871) . Irrespective of 
its exact form, (A J 2 ) describes a random-walk process, and 
is therefore closely related to the relaxation time. In antici¬ 
pation of the eventual necessity of introducing such approxi¬ 
mations, we forgo from the outset the attempt to write down 
explicit expressions for the diffusion coefficients. Instead, we 

1 We use the notation N(x;y), where x stands for any argument or set of 
arguments (scalar of vector) and y is a parameter, to denote the stellar number 
density per dx interval at y. The units of N are the inverse of those of x. For 
example, N(e\ t) = f 0 m d JN(s, J;t) is the stellar number density per unit 
specific energy at time t. 
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use them to define the relaxation time t r as the time required 
to diffuse in J 2 by J^, 


t r 0 ) 


JM 

(AT 2 ) ’ 


(6) 


We then treat the relaxation time as a free parameter that char¬ 
acterizes the system’s typical timescale for the evolution of 
the DF, and whose value can be estimated by Eqs. (12511261 
below. For simplicity, the relaxation time t r is assumed to be 
independent of angular momentum 2 , but can generally be a 
function of energy. 

At steady state, the stellar current T (e) is independent of 
J, 


m = i 

W 2 t r 


N(e, J ) dN(s, J) 


J 


dJ 


Solving this equation yields 

N{e, J) = ~ 


2 J~t r 




-JlnJ + CJ. 


(7) 


( 8 ) 


The integration constants C and T that are determined by 
the boundary conditions N(e,Ji c ) = 0 and N(£,J m ) = 
Ni so (s, J m ). The isotropic DF is separable in e and J, 


N iao (e, J)dedJ = 2 ^ 2 °^^ ded J. (9) 

■'m l e J 

Applying these boundary conditions 3 to equation {SJ, the DF 
is given by 


2 N iso {e)J In (J/J lc ) 
Jl(e) HJm/JicV 


( 10 ) 


and the stellar current into the MBH per energy interval is 


■T(e) = , . T “°/7 , . ■ (11) 

In (Jm/Jic)t r {e) 

Note that the capture rate in the diffusive regime depends only 
logarithmically on the size of the loss cone. 

The prompt infall rate Y' p in the diffusive regime is then 
given by 


2.2. Keplerian orbits in a power-law cusp 

The MBH dominates the stellar potential within the radius 
of influence. 


Th 



(13) 


where a 2 is the ID stellar velocity dispersion far from the 
MBH. The mass enclosed within r/,. is roughly equal to M,. 
Various formation scenarios predict that the spatial stellar 
number density at r < ru should be approx imately a power 
law (e.g. Bahcall & Wolf ll975 Young?1980) 


n*(r) 


(3/2 - p)N h (r_y 3/2 ~ P 
4tt rl \r h J 


(14) 


where N h is the number of stars inside 77 ,. This corresponds 
to an energy distribution N(e)de ex e p_5 ( 2 de. A stellar cusp 
with 1 0 has been observed in the Galactic Center (Alexan¬ 
der 1929|; Genzel et al. 12003 ). Fora single mass population, it 
was shown by analytical considerations that p =1/4 (Bahcall 
& Wolf l 19761) . This has been confirmed recently by N-body 
simulations (Baumgardt et al. 12004a : Preto et al. 120041) . 

Mass segregation drives the heavy stars in the population 
to the center. The radial distribution of the different mass 
components can be then approximated as average power-laws, 
steeper (p > 0 ) for the heavier masses and flatter (p < 0 ) for 
the lower masses (Bahcall & Wolf Il97~7t Baumgardt et al. 
t2QQ4b). 

Typically, the diffusive regime is within the radius of influ¬ 
ence. We therefore assume from this point on that the stars 
move on Keplerian orbits (vj = GM m /r) in a power-law den¬ 
sity cusp. The stellar orbits are characterized by a semi-major 
axis a, eccentricity e, periapse r p and period P, 


a = 


GM. 

2e 


r p =a{ 1 - e), 


! = 1- 


J 2 

GM.a ’ 


P = 


2 ™ 3 / 2 

VGW.' 


(15) 


During most of the inspiral 1 — e <C 1, and the periapse can 
be approximated by r p ss J 2 /2CM This remains valid until 
the last phases of the inspiral. 

The prompt infall rate (Eq. O can be expressed in terms 
of the maximal semi-major axis. 


r p = 


de7V iso (e) 
l n (T TO / Jic)t r {e) 


( 12 ) 


T P = 


daNi so (a) 

1 Jlc)t r ((l) 


(16) 


where the energy e p separates the diffusive and full loss-cone 
regimes. A star samples all angular momenta Ji c < J < J m 
in a relaxation time, and it is promptly captured once J <Ji c . 
The total rate is therefore of order T p ~ Ni so (< a p )/t r (a p ), 
where a p is the typical radius associated with orbits of en ergy 
e p (a. p is the semi-major axis for Keplerian orbits, see 12.2> 
and Ni so (< a p ) is the number of stars within a p . The rate is 
logarithmically suppressed because of the diluted occupation 
of phase space near the loss cone. 

2 In general, if the relaxation time depends on r, this will introduce some 
dependence of t r on J. We do not consider this dependence here. 

3 Adiabatic MBH growth may lead to some anisotropy (Quinlan, Hemquist 
& Sigurdsson l 1995 1). Here we assume that far from the loss-cone the DF will 
be isotropic 


We will assume Keple rian orbits throughout most of this 
paper except for section (14. Ill , where we employ the general 
relativistic potential of the MBH. 

2.3. Slow inspiral 

The derivations of the conditions necessary for slow inspi¬ 
ral and of the inspiral rate follow closely those of prompt in¬ 
fall, but with two important differences. (1) The time to com¬ 
plete the inspiral is not the infall time I\ s , but rather to 3> -Po- 
(2) There is no contribution from the full loss-cone regime, 
where stars are scattered multiple times each orbit. This is be¬ 
cause inspiral in this regime would require that the very same 
star that was initially deflected into an eccentric orbit, be re¬ 
scattered back into it multiple times. The probability for this 


































4 


HOPMAN & ALEXANDER 


happening is effectively zero 4 . 

In analogy to the radial scale a p of prompt infall, which 
delimits the volume where stars can avoid scattering for a time 
Po < tj, and thus maintain their infall orbit until they reach 
the MBH, the inspiral criterion / (l < tj defines a critical radius 
a c <C (ip (or equivalently, a critical energy e c 3> e p ). Stars 
starting the inspiral from orbits with ao < a c (eo > £ c ) will 
complete it with high probability, whereas stars starting with 
ao > a c (e o < £ c ), will sample all J values before they spiral 
in significantly, regardless of Jo and ultimately either (1) fall 
in the MBH, (2) diffuse in energy to much wider orbits or (3) 
into the much tighter orbits of the diffusive regime. Since we 
assume a steady state DF, outcomes (2) and (3) represent a 
trivial, DF-preserving redistribution of stars in phase space, 
which does not affect the statistical properties of the system. 
Whether stars spiral in or fall in depends, statisticall y, on ly on 
ao (£o). We use below Monte Carlo simulations ( fl4.ll Fig. 
0 to estimate the inspiral probability function, S(ao), which 
describes the probability of completing inspiral when starting 
from an orbit with semi-major axis ao (S —>> 1 for ao <C a c , 
S —> 0 for ao a c ). The inspiral rate for stars of type s with 
number fraction f s is then 


r d aN(a)S(a) d aJV(a) 

ls Jo ln(J m /J Jc )f r (a) “ s Jo ln(J m /J ic )f r (a) ’ 

(17) 

where roughly S ( a c ) ~ 0.5. 

3. PARAMETER DEPENDENCE OF THE INSPIRAL RATE 

In this section we derive some analytical results for the in¬ 
spiral rate. In order to keep the arguments transparent, we 
neglect relativistic deviations from K epler ian motion. Rela¬ 
tivistic orbits are discussed in section d4. 1 k 
Consider a star of mass M* orbiting a MBH of mass M, on 
a bound Keplerian orbit with semi-major axis a and angular 
momentum J. When the star arrives at periapse, it loses some 
orbital energy A E by GW emission. As a result, the orbit 
shrinks and its energy increases. For highly eccentric orbits 
the periapse of the star is approximately constant during in¬ 
spiral in absence of scattering. We define the inspiral time to 
as the time it takes the initial energy £q to grow formally to 
infinity. If the energy loss per orbit is constant, then for e —> 1 


to = 


d£ 


£ o 


(d£/df) 


1 

A E 


or 


to{r p ,a) = 

For GW, A E is given by (Peters [1964) 


nOO 

/ d eP(e) = 

J co 

2e 0 Po 
AE ’ 

(18) 

2-k\J GM,a 


(19) 

AE 



AE 0W = 

5^2 M. 


-7/2 


where 


1 + — e 2 + — e 4 
/(e) = 24 96 


(1 + e) 7 / 2 


( 20 ) 


( 21 ) 


and rs = 2GM./c 2 is the Schwarzschild radius. 


4 This is to be contrasted with prompt infall from the full loss-cone regime, 
where a star can reach the MBH by being scattered once into the loss-cone 
just before crossing a p toward the MBH. 


During all but the last stages of the inspiral e ~ 1, in which 

case r p /r s = 4(J/J ic ) 2 and 


AEqw = E\{J/Ji c ) 


-7 


Ex = 


85n M+c 2 


3x 2 13 M. 
Gravity waves also carry angular momentum. 


. 16tt . 

A ./gw =-—5(e)- 

5 c 


5(e) = 


1 


le 2 


(1 + e) 2 ' 


( 22 ) 


(23) 

(24) 


Generally, the change in J in the course of inspiral is dom¬ 
inated by two-body scattering, and A Jgw can be neglected 
until a becomes very small. 

It is convenient to refer the timescales in the system to the 
relaxation time at the MBH radius of influence. 


, (M.V P(r h ) 
h p {mJ A/logA! 


(25) 


where Ai = M./M* (rs/r^,) 1 / 4 (Miralda-Escude & Gould 
120001) . and A p ~ 0.2 for p = 0 (Alexander & Hopman 12003 1. 
The relaxation time at any radius is then 


t r (a) = th (jr'j > (26) 

where we associate the typical relaxation time on an orbit with 
that at its semi-major axis. This is a good.approximation, 
since theoretical arguments (Bahcall & Wolf 1976 : 1977}) and 
simula tions (Freitag & Benz l2002t Baumgardt et al. [2004a 
2004 h: Preto et al. 120041) indicate that 0 < p < 0.25, and so 
t r is roughly independent of radius. The angular momentum 
relaxation time is 


tj 




(27) 


Dissipational inspiral takes place in the presence of two- 
body scattering. When to ~ tj, both effects have to be taken 
into account. It is useful to parametrize the relative impor¬ 
tance of dissipation and scattering by the dimensionless quan¬ 
tity 


s( J, a) 


to(J,a) 
tj(J,a ) 



(28) 


where we introduce the (^-independent) length scale 

, (8 VGW.Etfh 

dc = - 2 - 

which is of the same order as a c (Eq. 1301) and is <C r^. We 
define some critical value ,s cr it •O1 such that the inspiral is so 
rapid that the orbit is effectively decoupled from the perturba¬ 
tions. 

The three phases of inspiral can be classified by the value 
of s. In the “scattering phase” the star is far from the region 
in phase space where GW emission is efficient and s /g> 1. 
With time it may scatter to a lower-J orbit, enter the “tran¬ 
sition phase”, where s ~ 1, and start to spiral in. If it is not 
scattered into the MBH or to a wide orbit, it will eventually 
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reach the stage where s < s cr it- It will then enter the “dissipa¬ 
tion phase” where it spirals-in deterministically according to 
equations (1204424k Note that eventually the e—> 1 approxima¬ 
tion is no longer valid. 

Here we are mainly interested in understanding how the in¬ 
terplay between two-body scattering and energy dissipation in 
the first two phases sets the initial conditions for the GW emis¬ 
sion in the final phase. It should be emphasized that the onset 
of the dissipation phase does not necessarily coincide with the 
emission of detectable GW. For example, while a star is well 
into the dissipation phase by the time s 5 cr it ^ 10 , the 

orbit has still to decay substantially before the GW frequency 
becomes high enough to be detected by LISA. 

We derive an analytical order of magnitude estimate for the 
critical semi-major axis a c by associating it with s = 1 orbits 
in the transition phase. Since s falls steeply with J, we set 
J = Jic and solve s(Ji c , a c ) = 1 for a c , obtaining 


increases the prompt disruption rate, but will not enhance the 
rate of inspiral events. 

The dependence of the GW event rate on the mass of the 
MBH can be estimated from Eq. ED and the empirical Alt¬ 
er relation 

M. = 1.3X10 »M 0 ( 5sn l- FI ) a . (34) 

where /3~4 (Ferrarese & Merritt f2000 : Gebhardt et al. 200 0:: 
Tremaine et al. 20021 Note that the M,-a relation im¬ 
plies that the stellar number density at the radius of influ¬ 
ence, rih is larger for lighter MBHs: for example, for (3 = 4, 

rih oc Nhr^ 3 cx M, 1 ^ 2 , where we assumed that Nh oc M,; the 
consequences of this for the dependence of the rate for prompt 
tidal disruptions on the MBH mass were discussed by Wang 
& Merritt (20044. 

The GW event rate depends on Al, as 


etc _ / d c \ 

~r h ) 


Th 


3/(3 —2 p) 


(30) 


The MC simulations below ( 49. 1 1 confirm that this analyti¬ 
cal estimate corresponds within a factor of order unity to the 
semi-major axis where the inspiral probability S(a c ) ~ 0.5, 
for a wide range of masses (see table El). Expression (1171 for 
the rate can then be approximated by 


fsN h /Oc\ 

Ih l n ['Im(etc)/'Izc] yt'h J 


3 / 2-2 p 


(31) 


where Nh is the number of stars within rh- 
The dependence of the inspiral rate on p (at fixed d c /rh and 
neglecting the logarithmic terms) can be examined by writing 

r i^g(p)f s {N h /t h ), where 


Ti oc M. 3//3 1 . (35) 

This dependence is weak, e.g. T* oc M, 1 ^ 4 for /3 = 4. Thus, 
the rate becomes higher for lower mass MBHs. If Al, ~ 
10 3 M 0 IBHs indeed exist in stellar clusters and the M,-a 
relation can be extrapolated to these masses, then they may 
be more likely to capture stars than MBHs. This, however, 
does not necessarily translate into more GW sources. The 
strain h of GW decreases with the mass of the MBH, so that 
these sources have to be closer by in order to observe their 
GW emission. Another restriction is that for IBHs the tidal 
force is so strong that white dwarfs are tidally disrupted well 
outside the event horizon, which precludes them from being 
LISA sources. These issues are further discussed in 0 

4. ORBITAL EVOLUTION WITH DISSIPATION AND SCATTERING 


/ j \ 3(3—4p)/2(3 —2p) 

*> = (£) ■ 02) 
The pre-factor q(v) grows with p over the relevant range (see 
also Ivanov i2002 ). This reflects the fact that the inspiral rate 
is determined by the number of stars within a c , rather than 
the total number of stars within rh . The concentration of 
the cusp increases with p, so that there are more stars within 
a c . This result suggests qualitatively that in a mass-segregated 
population, the heavier stars (higher p) will have an enhanced 
GW event rate compared to the light stars (lower p). 

(■J-'rom equations (12914321 it follows that 

Ti oc f^ 2p/(3 “ 2p) . (33) 

Since 0 for typical stellar cusps around MBHs, this means 
that the inspiral rate is nearly independent of the relaxation 
time for such cusps. This counter-intuitive result reflects the 
near balance between two competing effects. When scatter¬ 
ing is more efficient, stars are supplied to inspiral orbits at 
a higher rate, but are also scattered off them prematurely at 
a higher rate, so the volume of the diffusive regime, which 
contributes stars to the inspiral (~ a 3 ), decreases. This is in 
contrast to the prompt disruption rate T p , which increases as 
the relaxation time becomes shorter 5 . It then follows that en¬ 
hanced scattering, such as by massive perturbers (e.g. clus¬ 
ters, giant molecular clouds; Zhao, Haehnelt & Rees 12001 

5 The prompt disruption rate is F p ~ N{< a p ) {rs / a p ) / P(a v ) (Syer & 
Ulmer IT99S Eq. 10). For prompt disruption a p ~ t^ 5-2p > (Alexander & 
Hopman|2003j), and so the rate scales as T p ^i“( 2 + 2j 9/L -2 P) 


We present three different methods for analyzing inspiral 
in the presence of scattering. The first approach is based on 
Monte Carlo (MC) simulations, which follow a star on a rel¬ 
ativistic orbit, described by e and J, and add small pertur¬ 
bations to simulate energy dissipation and random two-body 
scattering. The second approach consists of direct numerical 
integration of the time dependent diffusion-dissipation equa¬ 
tion. The third approach is a heuristic semi-analytical effec¬ 
tive model that can describe the “effective” trajectory of a star 
through phase space, as well as the statistical properties of an 
ensemble of such trajectories. 

The three approaches are complementary. The MC simu¬ 
lations allow a direct realization of the micro-physics of the 
system, since they follow the perturbed orbits of individual 
stars. They also offer much flexibility in setting the initial 
conditions of the numerical experiments, but the results are 
subject to statistical noise and are not easy to generalize. The 
diffusion-dissipation equation on the other hand deals directly 
with the DF, and allows an analytical formulation of the prob¬ 
lem in terms of partial differential equations, which are solved 
numerically. However, computational limitations do not allow 
covering as large a dynamical range as in the MC simulations. 
Finally, the heuristic effective model has the advantage of its 
intuitive directness and relative simplicity of use. We find 
that all three methods give the same results for the same un¬ 
derlying assumptions (Fig. [2}- This inspires confidence in the 
robustness of the analysis. 

To compare the three methods, we stop the simulation at 
the point where s = s cr it = 10~ 3 , and we plot the DF of the 
angular momenta of the stars. From that point on the stars 













6 


HOPMAN & ALEXANDER 


are effectively decoupled from the cluster and can spiral in 
undisturbed. We then use the MC method, which can be easily 
extended to follow the stars in the dissipation phase (section 
0, to find the DF of the eccentricities of stars which enter the 
LISA band. 


4.1. Monte Carlo simulations 

The MC simulations generally follow the scheme used by 
Hils & Bender (fi995h to study the event rate of GW. The star 
starts on an initial orbit with (eo, Jo) such that dissipation by 
GW emission is negligible. The initial value of J is not of 
importance as the angular momentum is quickly randomized 
in the first few steps of the simulation. Every orbital period 
P(e), the energy and angular momentum are modified by 8e 
and 5J, and the orbital period and periapse are recalculated 
(this diffusion approach is justified as long as P <C tj). The 
simulation stops when the star decays to an orbit with s < 
.s'crit =10 3 or when J < J; c and it falls in the MBH (escape to 
a less bound orbit is not an option here since energy r elax ation 
is neglected and only dissipation is considered; see ' 12.3k 

When stars reach high eccentricities as a result of scatter¬ 
ing, the periapse approaches the Schwarzschild radius to the 
point where the Newtonian approximation breaks down.The 
MC simulations take this into account by integrating the or¬ 
bits in the relativistic potential of a Schwarzschild MBH. The 
periapse of the star moving in a Schwarzschild spacetime is 
related to its angular momentum by one of the three roots of 
the equation 

=(<= 2 -^)(= 2 + ^)aV OR(r ,J). (36) 

The term on left hand side of this equation is the squared spe¬ 
cific relativistic energy of the star (including its rest mass), 
while the right hand side is the effective GR potential. A star 
on a bound orbit is not captured by the MBH as long as equa¬ 
tion d36j has three real roots for r. The smallest root is ir¬ 
relevant for our purposes. The intermediate root (the turning- 
point) is the periapse r p of the orbit, and the largest root is 
the apo-apse r tt . The semi-major axis a and eccentricity e are 
defined by a = (r p + r a )/2 and e = 1 — r p /a (see Eq. fl5l . 
For bound orbits the loss cone is a very weak function of 
energy and is well approximated by Eq. m 

For given values of E and J, the GR periapse is smaller 
than the Newtonian one, and therefore the orbits can be more 
eccentric, and the dissipation can be stronger than implied by 
the Newtonian approximation. The Keplerian relations be¬ 
tween energy, angular momentum and the orbital parameters 
(Basil 5> are replaced by 


E 


2 _ 

GR — 


J 2 = 


(q — 2 — 2 e)(q — 2 + 2e) 

(37) 

q(q - 3 - e 2 ) 

q 2 (GM,\ 2 

q — 3 — e 2 \ c J 

(38) 


where q = 2(1 — e 2 )a/rs (Cutler, Kennefick & Poisson 
11994 ), The condition for a non-plunging orbit can also be ex¬ 
pressed in terms of the eccentricity and the semi-major axis, 
2(1 - e 2 )a = (6 + 2e)rs (Cutler et al. 1 1994 . This corre¬ 
sponds to a maximal eccentricity for a star on a non-plunging 
orbit 

e m ax(a) = -tt + V (?’s/2a) 2 - 3 r s /a + 1, (39) 


which increases with a. If cll = (Pl/d^GM,) 1 ^ 3 is the 
maximal semi-major axis a star may have in order to be de¬ 
tectable by LISA , the maximal eccentricity of a star detectable 
by LISA is e max (a L ). 

The step in J-space per orbit is the sum of three terms 

8J(e, J) = A 1 J scat (e, J) + xA 2 ^scat (e, J) - AJgw(g J) • 

(40) 

The first a nd second terms represent two-body scattering 
(Eas. 1213161 . with 6 Ai J scat = (AJ) P= J'^Pfe) / {2t r J) and 
A 2 J scat = y/{A J 2 )P = sjP(s)/t r J m (e ) = y/P{e)/tjJ. 
The random variable x takes the values ±1 with equal prob¬ 
abilities. The third term is the deterministic angular momen¬ 
tum loss due to GW emission (Eq. El .The energy step per 
orbit is deterministic (diffusion in energy is neglected) 

6E{e,J) = AE GW (e,J). (41) 

In order to increase the speed of the simulation we use an 
adaptive time-step. We checked for some cases that taking a 
smaller time-step does not affect the results. 

The DF of inspiraling stars is generated by running many 
simulations (typically 3 x 10 4 ) of stellar trajectories through 
phase space with the same initial value for Jq for given initial 
semi-major axis oo, but with different random perturbations. 
We verified that the initial value of Jq is irrelevant, as long as 
Jo 3> Jic ■ We record the fraction of stars that avoid falling in 
the MBH, S(ao), and the value of J at s cr it for those stars that 
reach the dissipation-dominated phase, thereby obtaining the 
DF W(J\ a o). This is repeated for a range of a o values. The 
integrated DF over all cusp stars, W ( J ), is then obtained by 
taking the average of all DFs, weighted by N(ao)S(ao)dao 

(cfEq.O- 

4.2. Diffusion equation with GW dissipation 

The DF of the inspiraling stars at the onset of the dissipa¬ 
tion dominated phase can be obtained directly by solving the 
diffusion-dissipation equation. This approach was taken by 
Ivanov J2001 - who included a GW dissipation term and ob¬ 
tained analytic expressions for the GW event rate for J Jj c 
under various simplifying assumptions. Here we are inter¬ 
ested DF of stars very near the loss-cone, and so we inte¬ 
grate the diffusion-dissipation equation numerically with two 
simplifying assumptions. (1) We neglect the drift term in the 
diffusion-dissipation equation. This can be justified by noting 
that the drift grows linearly with time, 5J\{t) = J^ n t/{2t r J), 
while the change due to the random walk grows as jj 2 (i) = 
(t/tr) 1 ! 2 J m . For a star with initial angular momentum J the 
drift becomes important only for t > tdrift = 4(J/J m ) 2 f r = 
4 tj. Since we are interested in the timescales t ~ to < tj, 
the drift is only a small correction. (2) We assume Keplerian 
orbits. These approximations are validated by the very good 
agreement with the MC simulation results (section fOl . 

With these two assumptions, the diffusion-dissipation is (cf 
Eqs.00 


dN(e,J;t) 

dt 


Jm(e) d 2 

2 t r dJ 2 


N(e,J;t) 


6 The drift term Ai J SC at represents a bias to scatter away from the MBH 
due to the 2D character of the direction of the velocity vector v. This can 
be expressed geometrically by considering a circle of radius Au sca t centered 
on v (v + Av sca tis the change per orbit due to scattering). This circle is 
intersected by a second circle of radius v that passes through v and is centered 
on the radius vector to the MBH. The section of the first circle that leads away 
from the MBH is slightly larger than the section leading toward it. 
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start at some given initial semi-major axis ao and angular mo¬ 
mentum Jo, and follow the star during its first two stages of 
inspiral, i.e., until the moment that s < s cr it. 

The DF of stars which reach s = s cr it is determined by 
the evolution of the orbital parameters (J, a) in the region 
where energy dissipation is efficient. Because of the strong 
dependence of the energy A E which is dissipated per orbit 
on angular momentum (or, equivalently, periapse), this is a 
small region AJ in J-space, the size of which is of order 
of the loss-cone, A J ~ J/,.. The time At it typically takes 
a star with semi-major axis a o performing a random walk 
in J to cross this region is Af = [AJ/J m (ao)] 2 th ( t r = th 
assumed). This can be used to introduce an effective J- 
“velocity” AJ/At = -J/ n (a 0 )/Ji c t h . 

The effective velocity of the semi-major axis is given by 


Fig. 1.— Examples of tracks of stars from the MC simulation (solid lines), 
and a track from the effective model (dashed line) with £ = 1. The effective 
track represents the stochastic tracks well. The tracks extend to the point 
where s = 10 —3 . The relevant parameters are given in the first line of table 
El . i.e. inspiral of a white dwarf into a MBH. 


+ 


d_ 

de 


EGw{J)N(e,J-,t) , 


(42) 


da( J, a) _ da • 2a 2 AE G w(J) 

d t ~ dE ~ GM . P(a) ' 

These two equation define a deterministic time evolution in 
(J, a) space from given initial values (J(0) = Jo, a(0) = ao), 
to a final point ( Jf, a/), where s = s cr it- To recover W( J; ao) 
at s cr it, we introduce some scatter in the effective velocities 
dJ/dt, 


where Eqw = A /few /P is the rate at which energy is lost to 
GW. The first term accounts for diffusion in J-space and the 
second represents the energy dissipated by GW emission. As 
with the MC simulations, the diffusion coefficient (AJ 2 ) = 
Jm/t r is an input parameter rather than resulting from a self- 
consistent calculation. 

The initial conditions consist of an isotropic cusp DF, 
N(s, J) oc e v ~ 5 / 2 J and the boundary condition that the DF 
vanish on the loci s = s cr it = 10 -3 and J = J; c in the (e, J)- 
grid. The initial DF is evolved in time over the (e, J)-grid 
until t ~ tj, when a relaxed steady state is achieved. The inte¬ 
gration is done using a Forward Time Centered Space (FTCS) 
representation of the diffusion term, and an ’upwind’ differ¬ 
ential scheme for the dissipation term (see e.g. Press et al. 
um. After each time step (which are chosen small enough 
to obey the Courant condition) the DF is re-normalized so that 
captured stars are replenished. 

After relaxation, the DF is non-zero up to the boundary 
s = s cr it, as stars are redistributed in phase space by diffu¬ 
sion at large J and by dissipation at small J. The DF at s cr it is 
then extracted to construct the angular momentum DF, W(J). 
Although the stars continue their trajectory in phase space be¬ 
yond s cr it all the way to the last stable orbit, this does not 
change W(J) because the rapid energy dissipation does not 
allow much J-diffusion. 


<±1 = 

dt th 


(44) 


where j denotes angular momentum in terms of J/,,, y is 
drawn from the positive wing of the normalized Gaussian dis¬ 
tribution, Gi(y), and £ ~ 1 (a free parameter) is the width 
of the distribution, whi ch can vary depending on the system’s 
parameters. Equations 6H can be solved analytically. 


th 


a(t) = 


(J e /a 0 ) 3/2 l 2 

6 Cyj 6 


ao . 


(45) 

The stopping condition s(j, a) = s cr j t gives an additional 
relation between j and a, so that the final values are related 
by af = s CT njf d c . This ties the initial and final values 

through the equations of motions 


y(jf) 


(6£) 1 {d c /a 0 ) 3 / 2 j f 6 

1 - s crit(^c/ a o) 1 // 2 J / 5/3 


(46) 


The relation between the initial distribution of j- 
velocities and the final j distribution at is W (j /; «o) = 
G i (y) [d y(jf)/ d j /] . The integrated DF over all the cusp stars 
is 


4.3. Analytical model of effective orbit 

Stars follow complicated stochastic tracks in (J, a) space; 
due to scattering they move back and forth in J, while they 
always drift to smaller a due to GW dissipation (energy dif¬ 
fusion by scattering is neglected). The drift rate depends 
strongly on J. Figure ffl shows some examples of stellar 
tracks in phase space, taken from the MC simulations. 

In this section we propose a heuristic analytical model, 
which captures some essential results of the MC simulations, 
while providing a more intuitive understanding. 

In this model we define the equations of motion of an “ef¬ 
fective track” of a star, which is deterministic and can be 
solved analytically. With this method we follow stars that 


W(jf)dj f = 


f ae N(a)W(jf, a) da 
f ac N(a) da 


(47) 


4.4. Comparison of the different methods 

We compare the results of the MC simulation, the diffusion- 
dissipation equation and an integration of equations m and 
63. In all cases the calculation is stopped when s = 10 3 ; 
at that point the eccentricity is still approximately unity, but 
scattering becomes entirely negligible so that even in the MC 
simulation the quantitie s evolve essentially deterministically 
according to equations J20I24> . The calculations assume a 
cusp with slope p = 0, M . = 3x 1O 6 M 0 and M* = 0.6 
(corresponding to WDs). 
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Fig. 2.— A comparison of the DFs of angular momenta for inspiraling stars 
in a p = 0 stellar cusp, derived from (1) MC si mulations, (2) direct integration 
of t he diffusion-dissipation equation (Ea. 1421 and (3) the effective model (Eq. 
EH with £ = 0.6). The DFs are normalized over the displayed J-range. The 
DFs are shown at s = 10“ 3 , when the dynamics are dominated by dissipation 
and scattering no longer plays a role, but the orbits are still well outside the 
LISA band, P Pl and e —► 1. This DF sets the initial conditions for 
the dissipation phase. The DF of the isotropic parent population (no sink, no 
dissipation; Eq.[9j is also shown for contrast. Unlike the DF of the inspiraling 
stars, it is dominated by high- J stars. 


Table 1 

Comparison of assumptions in the 3 methods 


Method 

GR potential Ai J scat 

A Jqw Stop at s crit 

Monte Carlo 

yes yes 

yes 

no“ 

Diffusion/Dissipation 

no no 

no 

yes 

Effective track 

no no 

no 

yes 

a For the purpose of comparison to the other two methods, 
the MC simulation was terminated at s C rit- 



Figure o shows the very good correspondence between the 
three approaches, whose underlying assumptions and approx¬ 
imations are summarized in table m. One important con¬ 
clusion is that all stars enter the dissipation phase with very 
small angular momenta, 1 < J/ J; c < 2- For a suitable choice 
of £ ~ 1, the actual complicated stochastic tracks are well 
mimicked by the tracks of the “effective stars”. The semi- 
analytical approach not only identifies the correct scale of an¬ 
gular momentum of inspiraling stars, but reproduces the DF. 
Its practical value lies in the relative ease of its application 
compared to the time consuming MC simulations or the inte¬ 
gration of the diffusion-dissipation equation. 

5. INSPIRAL RATES AND DISTRIBUTION OF ORBITAL 
PARAMETERS 

We now proceed to apply the MC simulation to different 
types of COs inspiraling into a 3x 10 6 M 0 MBH in a galactic 
center, and into a 10 3 M 0 IBH in a stellar cluster. We find 
from the simulations the critical semi-major axis a c and cal¬ 
culate the DF of the eccentricities of LISA sources. We stop 
the simulations when the orbital period falls below the longest 
period detectable by LISA , Pl = 10 4 s, which corresponds to 
a semi-major axis of 

a L = 23.5 M 6 - 2/3 r s , (48) 


Table 2 

Parameters for stellar populations and inspiral 
RATES 


Star 

M* 

th 

rh 

V 

fs 

a c 


F 


[Mq] 

[Gyr] 

[pc] 



[mpc] 

[mpc] 

[Gyr" 1 ] 

WD“ 

0.6 

10 

2 

0 

0.1 

30 

20 

7.8 

NS“ 

1.4 

5 

2 

0 

0.01 

40 

20 

1.8 

BH a 

10 

1 

2 

1/4 1(—3) 

20 

10 

4.7 

NS 6 

1.4 

1 (—3) 0.05 

0 

0.01 

2 

1 

4.3 

BH b 

10 

1(—3) 0.05 

1/4 

K-3) 

5 

2 

6.7 

° MBH with M. = 3x 10' 

3 Mq 





6 IBH with M. = 10 3 

Mq. 






c From equation <301 
d From the MC simulations 







where M, = 10 6 Mq M 0 .We record the eccentricity at that 
point and construct the DF of the eccentricity at Pl- It is 
straightforward to integrate the orbits in the eccentricity his¬ 
tograms backward (forward) in time to larger (smaller) values 
of P, see e.g. Barack & Cutler ( 120031) . 

The stars are distributed according to a powerlaw distribu¬ 
tion with different values for p. The total number of stars 
within rh was assumed to be Nh = 2 M,/Mq, with different 
number fractions for the respective species. See table □ for 
the assumed parameters of the stellar populations. 

5.1. Massive black holes in galactic centers 

We assume a M, = 3 x 10 6 M 0 MBH as a representative 
example and make the simplifying assumption that the MBH 
is not spinning. Real MBHs probably have non-zero angu¬ 
lar momentum (e.g. see for evidence of spin of the MBH in 
our Galactic Center Genzel et al. l2003f) . An important qual¬ 
itative difference is that in the Schwarzschild case the eccen¬ 
tricity of a GW-emitting star decreases monotonically with 
time, which this is not so for non-zero angular momentum 
(e.g. Glampedakis et al. {20021 . We conducted several MC 
simulations with Kerr metric orbits and verified that in spite 
of changes in details, our overall results hold. 

The tidal field of MBHs in this mass range disrupts main se¬ 
quence stars before they can become detectable LISA sources. 
A possible exception is our own Galactic Center, where the 
weak GW emission from very low mass main sequence stars 
(which are the densest and so the most robust against tidal dis¬ 
ruption) may be detected because of their proximity (Freitag 
12003 ). However, here we only consider GW from COs. Ta¬ 
ble 0 summarizes the parameters assumed for the properties 
of white dwarfs (WDs), neutron stars (NSs) and stellar mass 
black holes (BHs). The BHs are assumed to be strongly seg¬ 
regated in a steep cusp because of their much heavier mass 
(Bahcall & Wo If 119771) . The total (dark) mass in COs near 
MBHs is not known, but in future it may be constrained by de¬ 
viations from Keplerian motions of luminous stars very near 
the MBH in the Galactic center (Mouawad et al. l2004h . 

Fig. m shows the normalized inspiral probability func¬ 
tion S s (a o) , where s stands for WD, NS, and BH. The crit¬ 
ical semi-major axis is similar for WDs and NSs, but much 
smaller for BHs, because of the smaller relaxation time and 
the higher central concentration due to mass segregation. The 
functions S s (a) are used to calculate the inspiral rates (Eq. 
113 that are listed in table The rate for WDs is highest, 
Fwd = 7.8 Gyr -1 . We find that because of mass segrega- 
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Fig. 3.— Dependence of the normalized fraction of inspiral events on initial 
semi-major axis, for the case of a MBH of M* = 3 x 10 6 Mq (asterisks) 
and an IBH (triangles). The solid line is for white dwarfs, the dashed lines 
for neutron stars, and the dotted lines for stellar mass black holes. For the 
parameters of these cases see table 


tion, the rate for stellar BHs is of the same order of magnitude 
as that for WDs, Tbh = 4.7 Gyr -1 although their number 
fraction is lower by two orders of magnitude. The hierarchy 
of rates in table 0 is not, of course, necessarily the same as 
cosmic rates that LISA will observe; NSs and BHs are more 
massive than WDs, and can be observed at larger distances. 

Fig. 0 shows the eccentricity DFs of stars on orbits with 
P = Pl\ note that since a = cll is fixed, the orbits are fully 
determined by e. The DFs show a strong bias to large eccen¬ 
tricities. The maximal eccentricity possible for a star orbiting 
a 3x 10 6 Mq MBH in the LISA band is e m ax = 0.81 (EqJ |39}. 

It should be emphasized that the histogram in figure 0 can 
be obtained deterministically from the DF W (J; s cr it = 10~ 3 ) 
in figure 0, because the stars have already reached the point 
where scattering is negligible. This would not be the case had 
the MC simulation been terminated at .s cr ;t = 1 (e.g. Fre- 
itag 12003 1. since then significant scatterings would continue 
to redistribute the orbital parameters. We find however that 
the final DF (at P = Pl), which is obtained by integrating 
W (,/; s cr it = 1) forwa rd in ti me according to the GW dissi¬ 
pation equations (Eqs. I20I23> without scattering, is not much 
different from that shown in Fig. 0 . This coincidence is 
due to the fact that the stars with the largest eccentricities, 
which typically drop into the MBH, are replenished by the 
stars with slightly lower eccentricities. The main consequence 
of choosing s cr it too large is an overestimate of the total event 
rate; stars which actually fall into the MBH are erroneously 
counted as contributing to the GW event rate. Incidentally, 
even though stars that do not complete the inspiral are not in¬ 
dividually resolvable, they will contribute to the background 
noise in the LISA band (Barack & Cutlcr f2004 ). 

A premature neglect of the effects of scattering in previous 
studies (e.g. Freitag l200 1‘) probably explains in part why our 
derived rates are significantly lower. Those studies usually 
assumed that stars will spiral in without further perturbations 
once s = 1. We ran a simulation that was stopped at s cr it = 1 
instead of s cr j t = 1 0“ 3 . The event rate in that unrealistic case 
is about ~ 6 times higher. This does not explain all of the 
discrepancies, which are hard to track as different methods 
are used. The different stopping criterion may also explain 
why our rates are lower than those estimated by Hils & Bender 
( 1995 1. who used a method similar to ours, without specifying 



Fig. 4.— The probability DF of a CO entering the LISA band with eccen¬ 
tricity e for a MBH of mass M . = 3 X 10 6 Mq. The histograms represent 
WDs, NSs and BH (from broad to narrow). Note that the maximal possible 
eccentricity for which P < Pj, is emax = 0.81. See table 0 for the cusp 
parameters. 


the criterion for inspiral. 

5.2. Intermediate mass black holes in stellar clusters 

Unlike MBHs with masses M. > 10 6 Mq , there is little 
firm observational evidence at this time for the existence of 
IBHs with masses M, ~ 10 3 Mq (see review by Miller & 
Colbert 12004 1. However, there are plausible arguments argu¬ 
ing in favor of their existence. From a theoretical point of 
view, these objects are thought to form naturally, such as in 
population III remnants (Madau & Rees 1200 111 , or in a run¬ 
away merger of young jtars in dense stellar clusters (Portegies 
Zwart & McMillan l200ll Portegies Zwart et al. 12004 .1. From 
an observational perspective, IBHs may power some of the 
ultraluminous X-ray sources (e.g. Miller, Fabian, & Miller 
12001 - for example by tidal capture of a main sequence com¬ 
panion star (Hopman et al. 12004 1. 

For the purpose of estimating the orbital parameters of GW 
emitting stars, we assume that IBHs lie at the center of dense 
stellar clusters (see model parameters listed in table0. White 
dwarfs are disrupted by an IBH before entering the LISA band, 
so that only the most compact sources, neutron stars and stel¬ 
lar mass BHs can emit GW in the LISA frequency band. The 
same values for the stellar population fractions f s where taken 
here as for MBHs. We note that this is not necessarily the 
case. For example, N-body simulations indicate that a large 
fraction of BHs may be ejected in an early stage of the clus¬ 
ter’s life if a massive stellar object forms a tight binary with 
the IBH (Baumgardt et el. i2004b ii. 

Figure 0 shows the inspiral probability functions S s (a), 
and Fig. 0 shows the eccentricity histograms of stars on 
orbits with P = P^. The maximum eccentricity still observ¬ 
able by LISA is nearly unity, and the IBH case shows even 
more clearly the strong tendency towards large eccentricities. 
In general, this effect becomes more prominent for lighter 
BHs. The high eccentricity makes the GW signal highly non- 
monochromatic. The star spends most time at apo-apse, emit¬ 
ting a relatively weak, low frequency signal. At periapse 
short pulses of high frequency GW is emitted. For IBHs of 
M. ~ 10 3 Mq, the frequency v p of these short bursts at pe¬ 
riapse is of the order v p < 100 Hz. This is too high to be 
measurable by LISA , but may be measurable by ground based 
detectors such as LIGO or VIRGO. 
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0.986 0.988 0.99 0.992 0.994 0.996 0.998 


e 

Fig. 5.— The probability DF of a compact remnant entering the LISA band 
with eccentricity e for an IBH of mass M. = 10 3 Mq . Only NSs (broad) and 
BHs (narrow) are considered; WDs are probably disrupted by the tidal field. 
For an IBH, the maximal possible eccentricity for which P < Pl is nearly 
unity, emax = 0.998, all inspiraling stars are likely to have eccentricities 
close to the maximum value. See table 0 for the cusp parameters. 


As anticipated (Eq. E). the rate of inspiraling stars is com¬ 
parable to that of a MBH (Table [2}. However, due to the ex¬ 
tremely high eccentricities and small semi-major axes of the 
orbits, the GW emission is very efficient and they spiral in on 
a very short time scale (on the order of a year). 

6. SUMMARY AND DISCUSSION 

Stars spiraling into MBH due to the emission of GW are 
an important potential source for future GW detectors, such 
as LISA. The detection of the signal against the noise will 
be challenging, and requires pre-calculated wave-forms. The 
waveforms depend on the orbital parameters of the inspiraling 
stars. Our main goal in this study was to derive the distribu¬ 
tion of eccentricities of inspiraling stars. The orbital statis¬ 
tics of such stars reflect a competition between orbital decay 
through dissipation by the emission of GW, and scattering, 
which deflects stars into eccentric orbits but can also deflect 
them back to wider orbits or straight into the MBH. 

Inspiral is a slow process, and unless the stars start close 
enough to the MBH, they will be scattered off their orbit. We 
identified a critical length scale a c which demarcates the vol¬ 
ume from which inspiral is possible. We obtained an analyt¬ 
ical expression for the inspiral rate and showed that it is of 
the order of a few per Gyr per galaxy, much smaller than the 
rate for direct capture (Alexander & Hopman 1200 31. that it 
is nearly independent of the relaxation time for typical stellar 
cusps, and that it grows slowly with decreasing MBH mass 
(assuming the M. — a relation). Throughout we assumed 
a single powerlaw DF. Generalization to different profiles is 
straightforward. Qualitatively, the inspiral rate is determined 
by the stellar density near a c , and is not very sensitive to the 
exact profile far away at radii much smaller or much larger 
than a c . The rate of GW events depends on the number of 
COs inside a c (which is much smaller than the MBH radius 
of influence), and so mass-segregation can play an important 
role in enhancing the event rate by leading to a centrally con¬ 
centrated distribution of COs. 

We obtained a relatively low rate. One important reason 
for this is our stringent criterion for inspiral. This may also 
explain why our results deviate from Hils & Bender l 19951 . al¬ 
though they do not specify their precise criterion for inspiral. 
Comparison with other works are more complicated. Possi¬ 


ble differences may stem from different normalizations of the 
central density, different CO fractions, different criteria for 
inspiral, and the way mass segregation is treated. An essen¬ 
tial step in the future will be to analyze inspiral processes by 
N-body simulations, which are feasible already for small sys¬ 
tems (Baumgardt et al. 2004a H2004ht Preto et al. 120041) . 

The detection rate depends on the inspiral rates, but also on 
the mass, relaxation time, the orbital parameters (especially 
the eccentricity) and the details of the detection algorithm 
(Barack & Cutler H2003 T Here we provide a simple recipe 
to estimate the number of detectable sources. We assume that 
the various dependencies above can be expressed by an effec¬ 
tive strain ft. 

The strain of GW resulting from a star of mass M* = 
mM 0 orbiting a MBH of mass M. = 1Q 6 M@Mq 3> M* 
at a distance d = Gpc da pc , on a circular orbit of orbital fre¬ 
quency v = 10 ~ 3 s _1 u_ 3, is given by 


ft = 1.7 x Hr 23 u 2 {j 3 rf g p c Mq ^ 3 m . (49) 

(e.g., Sigurdsson & Rees ll997l ). 

The cosmic rate of LISA events is given by 


r 


tot 



dM,-jjpT i(M,)V (M,), 
dM. 


(50) 


where dn./dM. is the number of MBHs per unit mass, per 
unit volume, and V (M.) is the volume in which stars can be 
observed by LISA. Aller & Richstone (12002 1) used the M. — 
a relation to “weigh” the MBHs. The spectrum is roughly 
approximated by 


dn, 

dM. 


10 ' 


( 1 


\10 6 Mq 


Mg 1 Gpc" 


(51) 


The LISA sensitivity curve at v = 10 -3 is ft ~ 10 -23 for 
one year of observation with signal to noise ratio S/N=l (see 
e.g. http://www.srl.caltech.edu/lisa I. We adopt this value as a 
representative detection threshold for ft. If the effective strain 
has to be be at least 10 -23 ft_23 for the star to be observable, 
then, for a Euclidean Universe, 


V(M.) = 20.6 Gpc 3 ftl^z/ 2 3 to 3 M 6 2 . (52) 

Finally, the rate per MBH is 

r i(M.) = r;(3 x 1O 6 M 0 )(3M 6 ) -1 / 4 , (53) 

where we used the expression for the dependence of the in¬ 
spiral rate on the MBH mass, equation <I35> . with f3 = 4. 
The expression can be calibrated with the MC results for a 
3 x 1O 6 M 0 MBH. 

Integrating from Mq = 0.5 to Mq = 5 gives 


r tot = 1-5 yr 


T,(3 x 1O 6 M 0 ) 


Gyr 


-l 


t 3 ft: 3 3 u 2 3 . 


(54) 


The number of sources which LISA can observe at any mo¬ 
ment is estimated by Aft = Ttot x fi(e, ol), where fi is the 
time a star with eccentricity e spends in the LISA band before 
being swallowed; here e is the average eccentricity of stars 
entering the LISA band. 

For example, our calculations for WD inspiral indicate that 
Ti(3 x 1O 6 M 0 ) = 7.8 Gyr -1 . For a circular orbit with a 
period of P = 10 3 s and Mq = 1, the inspiral time is Q, = 
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54 yr, in which case the number of detectable sources would 
be A~ 130/^23^—3- 

We would like to emphasize that this estimate is to be 
treated with caution. The number of detectable sources de¬ 
pends strongly on the assumptions. In particular, GWs from 
eccentric sources are not monochromatic and may be harder 
to detect. The analysis in this paper can be used for a more 
detailed analysis of the number of detectable LISA sources. 

We used three complementary methods to model the in¬ 
spiral process, MC simulations, numerical solution of the 
diffusion/dissipation equation and an analytic effective orbit 
model. We followed the evolution of the orbital properties of 
the inspiraling stars to the point where they decoupled from 
the scattering. All three methods were in excellent agree¬ 
ment. We find that the distribution of orbital angular momenta 
is strongly peaked near J > Ji c , with the detailed form of 
the distribution depending somewhat on the parameters of the 
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